Spectral theory of metastability and extinction in a branching-annihilation reaction 
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We apply the spectral method, recently developed by the authors, to calculate the statistics of 
a reaction-limited multi-step birth-death process, or chemical reaction, that includes as elementary 
steps branching A — > 2A and annihilation 2A — > 0. The spectral method employs the generating 
function technique in conjunction with the Sturm-Liouville theory of linear differential operators. 
We focus on the limit when the branching rate is much higher than the annihilation rate, and obtain 
accurate analytical results for the complete probability distribution (including large deviations) of 
the metastable long-lived state, and for the extinction time statistics. The analytical results are in 
very good agreement with numerical calculations. Furthermore, we use this example to settle the 
issue of the "lacking" boundary condition in the spectral formulation. 
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PACS numbers: 05.40.-a, 02.50.Ey, 87.23.Cc, 82.20.- 



I. INTRODUCTION 



Statistics of rare events, or large deviations, in chem- 
ical reactions and systems of birth-death type have at- 
tracted a great deal of interest in many areas of sci- 
ence including physics, chemistry, astrochemistry, epi- 
demiology, population biology, cell biochemistry, etc. 

EE I SI, 1 ISllIiEiB Large 
deviations become of vital importance when discrete 
(non-continuum) nature of a population of "particles" 
(molecules, bacteria, cells, animals or even humans) 
drives it to extinction. A standard way of putting the 
discreteness of particles into theory is the master equa- 
tion 0,12] which describes the evolution of the probability 
of having a certain number of particles of each type at 
time t. The master equation is rarely soluble analyti- 
cally, and various approximations are in use [l|, ■ One 
widely used approximation is the Fokker-Planck equation 
which usually gives accurate results in the regions around 
the peaks of the probability distribution, but fails in its 
description of large deviations, that is the distribution 
tails [H, [ijj • Not much is known beyond the Fokker- 
Planck description. In some particular cases (especially, 
for single-step birth-death processes) complete statistics, 
including large deviations, were determined by applying 
various a ppr oximations directly to the pertinent master 
equation \ul [H, [H, CLl EE 51 M, SH, M, H • A differ- 
ent group of approaches employs the generating function 
formalism P, [|J Q , see below. Here the master equation 
is transformed into a linear partial differential equation 
(PDE) for the generating function, and this PDE is an- 
alyzed/solved by various techni ques such as the method 
of second quantization [24], l25l. l26j or the more recent 
time-dependent WKB approximation [27| • Recently, 
we combined the generating function technique with the 
Sturm-Liouville theory of linear differential operators and 
developed a spectral theory of rare events 28, 2jj. In this 
theory the problem of computing the complete statistics 
of (not necessarily single-step) birth-death systems re- 
duces to solving an eigenvalue problem for a linear differ- 
ential operator, the coefficients of which are determined 



by the reaction rates. 

In this paper we apply the spectral method to the 
paradigmatic problem of branching A + X — > 2X and 
annihilation X + X — > E, where A and E are fixed. This 
multi-step single-species birth-death process describes, 
for example, chemical oxidation reactions [HI, If 
the branching rate is much higher than the annihilation 
rate (the case we will be mostly interested in throughout 
the paper), a long-lived metastable, or quasi-stationary, 
state exists where the two processes (almost) balance 
each other. Still, this long-lived state slowly decays with 
time, because a sufficiently large fluctuation ultimately 
brings the system into the absorbing state of no parti- 
cles from which there is a zero probability of exiting. In 
this type of problems one is interested in the extinction 
time statistics and in the complete probability distribu- 
tion, including large deviations, of the quasi-stationary 
state (formally defined as the limiting distribution con- 
ditioned on non-extinction). Turner and Malek-Mansour 
[l8| calculated the mean extinction time in this system 
by solving a recursion equation for the extinction proba- 
bility. More recently, Elgart and Kamenev [l6[ revisited 
this problem in the light of their time-dependent WKB 
approximation for the generating function. Their insight- 
ful method readily yields an estimate of the mean extinc- 
tion time, but only up to a (significant) pre-exponential 
factor. The quasi-stationary distribution for this system 
has not been previously found, and calculating it will 
be our major objective. In the language of the spectral 
theory, the mean extinction time represents the inverse 
eigenvalue of the ground state, while the quasi-stationary 
distribution is derivable from the ground state eigenfunc- 
tion. 

The paradigmatic branching-annihilation problem, 
considered in this paper, has an additional value, as it 
helps settle one unresolved issue of the spectral theory. 
In the previous works [HI, [2i| we considered reactions 
that conserve parity of the particles. The parity conser- 
vation provides an additional boundary condition for the 
PDE for the generating function which ensures a closed 
formulation of the problem already at the stage of the 
time-dependent PDE. The branching-annihilation pro- 
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cess, considered in the present work, does not conserve 
parity. As we will show, the "lacking" boundary condi- 
tion emerges here (and in a host of other problems of this 
type) only at the stage of the Sturm-Liouville theory. 

Here is how we organize the rest of the paper. In 
Section II we apply the spectral method and reduce the 
governing master equations to a proper Sturm-Liouville 
problem. In Section III we employ a matched asymptotic 
expansion to approximately calculate the ground-state 
eigenvalue and cigcnfunction and obtain the long-time 
asymptotics of the generating function. This asymptotics 
is used in Section IV to extract the quasi-stationary prob- 
ability distribution and compare it with our numerical 
results. In Section V we calculate the mean extinction 
time and extinction probability distribution and compare 
these results with the previous work and with our numer- 
ics. Some final comments are presented in Section VI. 



II. GENERATING FUNCTION AND 
SPECTRAL FORMULATION 

We consider the branching and annihilation reactions 
A ^> 2A and 2A -% 0, where fi, A > are the rate con- 
stants. The (mean-field) rate equation for the number 
of particles n(t), dn/dt = Xn — fin 2 predicts a nontriv- 
ial attracting steady state n s = X/fi = Q. Fluctuations 
invalidate this mean-field result due to the existence of 
an absorbing state at n = 0. However, when O ^> 1, 
there exists a long-lived fluctuating metastable (or quasi- 
stationary) state, which slowly decays in time, implying 
a slow growth of the extinction probability. The statis- 
tics of this quasi-stationary state and of the extinction 
times are in the focus of our attention here. 

The master equation for the probability P n (t) to find 
n particles at time t can be written as 

J t P ^) = ^[(n + 2)(n+l)P n+2 (t)-n(n-l)P n (t)} 
+ X [(n - l)P n -i{t) ~ nP n (t)] , n > 1 , 



dt 



P (t) = fxP 2 (t). 



We introduce the generating function [l], 0, [1] 



G(x,t) 



5> n P„(t), 



(1) 



(2) 



where x is an auxiliary variable. Once G{x,t) is known, 
the probabilities P n (t) can be recovered from the Taylor 
expansion: 



Pn(t) = 



1 d n G{x,t) 



dx n 



(3) 



x=0 



By virtue of Eqs. ((2|) and ([3]), G(x, t) must be analytical, 
at all times, at x = 0. Equations |T]) and ^ yield a 
single PDE for G(x,t) [3]: 



dG 
~dt 



\(l-x 2 )^+Xx(x-l) 



dG 

dx 



(4) 



Conservation of probability yields one (universal) bound- 
ary condition for this parabolic PDE: G(l,t) = 1 [3lj ] - 
What is the second boundary condition? Note that 
G(x = — 1, t) must be bounded at all times, as it is equal 
to the difference between the sum of the probabilities to 
have an even number of particles and the sum of the prob- 
abilities to have an odd number of particles. Now, the 
steady-state solution of Eq. G s t{x) = G(x,t — > oo), 
that obeys the equation 



'-< I - ,r 2 ) d2Gst 

2 1 X > dx 2 



Ax (a 



1 



dGs_ 

dx 



= 0. 



(5) 



must also be bounded at x = — 1. Then Eq. ([5]) immedi- 
ately yields a second boundary condition: G' st (x)\ x= -i = 
0, where the prime stands for the x-derivative. Combined 
with G s t(l) = 1, this condition selects the steady state 
solution G s t{x) — 1 describing an empty state. 

Now let us introduce a new function g(x, t) — G(x, t) — 
G s t(x) = G(x,t) — 1 [which obeys Eq. dU) with a ho- 
mogenous boundary condition g(x = = and is 
bounded at x — —1], and look for separable solutions, 



— P -7fc* 



<Pk(x)- We obtain 



{l-x 2 )^{x)+2nx{x-l)ip' k (x) + 2E k tp k (x) =0, (6) 

where Ek = Ik/n-- One boundary condition is of course 
yfe(l) = 0. The second boundary condition comes from 
the demand that ifik(x) be bounded at x = —1. Then 
Eq. ^ yields a homogenous boundary condition 



2fV fc (-l) + £fcWc(-l) = 0, 



(7) 



for each k = 1,2,.... Notice that the eigenvalue Ek enters 
the boundary condition. Rewriting Eq. © in a self- 
adjoint form, 

[ip' k {x)e^(~2Qx)(l + x) 2n ]' + E k w(x)ip k (x) = 0, (8) 
with the weight function 



w(x) 



2e 



-2Qx 



(1 + x) 



2S2 



1 -X 2 



(9) 



we arrive at an eigenvalue problem of the Sturm-Liouville 
theory [3^ . Once the complete set of orthogonal eigen- 
functions ifk(x) and the respective real eigenvalues E kl 
k = 1,2,..., are calculated, one can write the exact so- 
lution of the time-dependent problem for G(x,t): 



G(x, t) = 1 + } a k (fk{x)e 



-fiE k t 



k=l 



where the amplitudes a k axe given by 

J_ AG(x, t = 0) — l]ip k (x)w(x)da 



ak = 



L 1 i fit(x)w(x)dx 



(10) 



(11) 



As all E k are positive, Eq. (fTOj) describes decay of ini- 
tially populated states k = 1,2,..., so the system ulti- 
mately approaches the empty state G{x,t — > oo) = 1. 
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Being mostly interested in the case of O 3> 1, we 
note that while the eigenvalues of the "excited states" 
E 2 ,E 3 ,... scale like 0(0) > 1 [s3|, the "ground state" 
eigenvalue E\ is exponentially small [Isj . Therefore, at 
sufficiently long times, /iQt = Xt 3> 1, the contribution 
from the excited states to G(x, t) becomes negligible, and 
we can write 



G(x, t) = 1 + a% <pi(x) e 



(12) 



So we need to calculate the ground-state eigenvalue E\, 
the eigenfunction (pi(x), and the amplitude a\. (Actu- 
ally, the eigenvalue E\ was calculated earlier [18j , but we 
will rederive it here.) Note that, as E\ is exponentially 
small, the boundary condition ([7]) for the ground state 
reduces, up to an exponentially small correction, to 



pi(-l)=0. 



(13) 



III. GROUND STATE CALCULATIONS 



Throughout the rest of the paper we assume O 3> 1. 
As .Ei is exponentially small, the last term in Eq. ^ is 
important only in a narrow boundary layer near x = 1, 
and we can solve Eq. ([6]) for ipi(x) = ip(x) by using a 
matched asymptotic expansion, see e.g. Ref. (34|. In 
the "bulk" region — 1 < x < 1 we can treat the last 
term in Eq. ^ perturbatively. In the zeroth order we 
put Ei = and arrive at the steady state equation (1 + 
x)ip"(x) — 2Clx<p'(x) = , whose (arbitrarily normalized) 
solution, bounded at x = —1, is ip^(x) — 1. Now we 
put ipb(x) = 1 + Scpb(x), where Sipb{x) -C 1, and obtain 
in the first order 



[6<p' b (x)e- 2Q *(l + x) 



\2m 



2Eie 



-2i7x 



-{l + x) 2n . (14) 



1 



Solving this equation, we obtain the bounded solution for 

ip b (x): 



Vb(x) = l-2Ei 



x e 2ns ds 



(1 + r) 2U e 



2Q„-2Cr 



-dr . 

(15) 

This solution, that obeys the boundary condition ([7]), is 
almost constant in the entire region — 1 < x < 1 except 
in the boundary layer near x = 1 (to be defined later on). 
To find the probabilities P n (t), we will need to calculate 
the derivatives of <fb(x) at x = 0. As long as 1 — x ^> I/O, 
we can neglect the r 2 term in the denominator of the 
inner integral in Eq. (|15p and obtain 



1 - 2Ei 



= 2fis 



ds 



o (1 + *) 



2!> 



Ei f e \ 2 ° 



(i 

i 

20s 



\2f2„-2n? 



dr 



(1 + s 

T[2fl + l,2Q(l + s)]}ds, 



|2<> 



{r[20 + i] 



(16) 



where T(a, z) = f °° s Q 1 e s (is is the incomplete 
Gamma function [3a|. Using the expansion [3^] 



r(a) - T(a,z) = ^2 



(17) 

we can evaluate the integral in Eq. (JT6J) and obtain 

E-\ ^r[2+j,-20]-r[2+j,-20(l+a;)] 



1 



\ " 



j!(20+i + l) 



(18) 



One can check that the perturbative solution in the bulk 
is valid [that is, 6<pb(x) <C 1] as long as 1 — x ^> I/O. 

In the boundary layer 1 - 1 < 1 we can disregard, in 
Eq. the (exponentially small) last term, and arrive at 
the same equation as before: (1 + x)(p"(x) — 2Ctxcp'(x) — 
. The solution obeying the required boundary condition 
at x = 1 is 



const x 



2Q.S 



'l + s)- 2n ds 



C 



1 



3 -2fi(l-ln2) 



(1+x) 



2Q 



(19) 



where C is a yet unknown constant. To find Ei and C we 
can match the asymptotes of the bulk and the boundary- 
layer solutions in the common region of their validity 
I/O <C 1 — x <C 1. Let us return to the first line of 
Eq. (|16p and evaluate f b (x) in this region. The inner 
integral receives the largest contribution from the vicinity 
of r = 0, while the outer integral receives the largest 
contribution from the vicinity of s = x. Therefore, we 
can extend the upper limit of the inner integral to infinity 
and obtain (see Appendix A) 



fb{x) 



1 - 2Ei 



„2<ix 



ds 



2i> 



1 



1 - 



2Ei^ 



X 



„2i1s 



'0 

2Ei^ • 



u 0- + S) 
2Slx 



(1 



ds 



r) 2n e- 2nr dr 



Q3/2 (! + a-jan ' 
Now, by matching Eqs. (Til?]) and (j2"0|) , we obtain 



Ei = 



OV2 
2^ 



,-an(i-ta2) and c = 1 



(20) 



(21) 



One can see that the ground state eigenvalue E\ is ex- 
ponentially small in 0. Equation (|2"Tj) yields the mean 
extinction time (/X-Ei) -1 (see Section V) which coincides 
with that obtained, by a different method, by Turner and 
Malek-Mansour [l|. 

Equations (fl~5| and (fT9|) yield the ground state eigen- 
function: 




for 1 - x > I/O , 
for 1 - x <C 1 . 



(22) 
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Now we use Eq. (|TT|) to calculate the amplitude a\ enter- 
ing Eq. (|12p . Let the initial number of particles be n , 



so G(x,t = 0) 



Evaluating the integrals, we notice 



that (i) the main contributions come from the bulk region 
1 — x ^S> I/O, and (ii) it suffices to take the eigenfunction 

ifib(x) in the zeroth order: tp^\x) ~ 1. Furthermore, 
when no 3> 1, the term x n ° under the integral in the nu- 
merator is negligible compared to 1. So, for n S> 1, the 
numerator and denominator are approximately equal to 
each other up to a minus sign. Therefore, a\ ~ — 1 (and 
independent of n Q ) which completes our solution (TT2]) for 
times fit ^> Q . 
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IV. STATISTICS OF THE QUASI-STATIONARY 
STATE AND ITS DECAY 

What is the average number of particles n(t) and the 
standard deviation a(t) at times fit ^> f2 _1 ? Using 
Eq. 42) and Eq. ([12]) with a x = -1, we obtain 



(23) 



n=0 



Furthermore, 



a 2 (t) = 



n 2 — n 2 



n=0 \n=0 / 

[d 2 xx G + d x G - (d x G) 2 ]\ x=1 

m 
2~ 



— + n 2 (i - e-^') 



,-fJ,E 1 t 



(24) 



where we have used for ipi (x) its boundary layer asymp- 
tote tfbi(x) (|19p with C = 1. At intermediate times 
Q^ 1 <C fit <C one obtains a weakly fluctuating quasi- 
stationary (metastable) state. Here the average number 
of particles, 



n. 



(25) 



coincides with the attracting point of the mean field the- 
ory, while the standard deviation, 



~~2 



1/2 



(26) 



coincides with that obtained from the Fokker-Planck 
approximation, see Appendix B. Note that a(t) from 
Eq. (|24p is a non-monotonic function. This stems from 
the fact that the quasi-stationary probability distribution 
around n ~ VL decays in time, whereas the extinction 
probability Po{t) grows. At times fiE\t <C 1 the stan- 
dard deviation a ~ ^/3Sl/2 corresponds to the unimodal 
quasi-stationary distribution around n ~ O, whereas at 
[iE]t ^> 1, a — > corresponds to the unimodal Kro- 
necker delta distribution at n — 0. At intermediate times, 
\xE\t ~ 1, the distribution is distinctly bi-modal. The 



FIG. 1: (Color online) (a) The average number of particles 
as a function of time at fit ^> I/O as described by Eq. (|23[) 
(the solid line) compared with the prediction from the rate 
equation n(t) ~ (the dashed line), for Q. — 30. Inset shows 
a blowup at intermediate times f2 _1 <C fit <C E^ 1 where 
the curves almost coincide, (b) The standard deviation from 
Eq. (|24p versus time for the same fl. 



maximum standard deviation <j max ~ 0/2 is obtained 
for e r^ Elt ~ 1/2. Figure ED shows the n[t) and a(t) de- 
pendences. 

Let us now proceed to calculating the complete prob- 
ability distribution P n (t) of the (slowly decaying) quasi- 
stationary state, conditional on non-extinction. For n — 
we obtain 



P (t) =G(x = 0,t) = l- e-» Elt 



(27) 



which, at fiE\t <C 1, is much less than unity. For n > 1 
Eqs. © and (0 yield 



P n (t) = - 



1 d n <p b (x) 



i! dx n 



(28) 



where ifb(x) should be taken from Eq. (|16p . After some 
algebra (see Appendix C), we obtain, for n > 1, 

PM= ^g|ga ■Bpn, w «,- M ).->«-, 

n 1 (2il + n) 

where \F\ (a, b, x) is the Kummer confluent hypergeomet- 
ric function [35| ■ To avoid excess of accuracy, we need to 
find the large O asymptotics of Eq. ([29]) . To that aim we 
use the identity [351 ] 



iFi(2f2,n+2f2, -20) 



T(n + 20) 

r(20)r(n) 



(l-s^ds (30) 



and consider separately two cases: n» 1 and n = 0(1). 

For n> 1, the integral in Eq. (f30|) can be evaluated by 
the saddle point method [34[. Denoting <I>(s) = —251s + 



5 



2fHn(s) + nln(l — s) we obtain 

2E1 V2^(2ny 



Pn(t) 



-i e 2n 



n\ ^20(1 

-2Q,\s 



s*) 2 + ns 2 

ln(s»)]+nln(l-s.)-/iSit ^g-^) 



where s„ = l+q — ((/ 2 + 2q) 1/ ' 2 is the solution of the saddle 
point equation $'(s) = 0, and q = nj (4S1). Equation (|3~lj) 
can be simplified in three limiting cases. In the high-n 
tail, n 3> O S> 1, we have s» ~ 257/n <C 1, and 



Pn(t) 



22f2-3/2 



20 

n 



n+2Q 



n-2U- 



In the low-n tail, 1 <C n <C O, we have s* 
and 



i-Vn7(2ny» 



>2f2-2 



20 

n 



ra/2+1/2 



%/2 — 2U—iiEit 



(33) 



Finally, for |n — 0| <C 0, s* 
we obtain 



1/2- (n- 0)/(60), and 



(34) 



For <C 1 this result describes a normal distribu- 

tion with mean O and variance 30/2, in agreement with 
Eqs. (|25|) and (|26[) and with the predictions from the 
Fokker-Planck equation, see Appendix B. 

Now we turn to the case of n = 0(1). Then it is always 
n <C O. Here it is convenient to rewrite the integral in 
Eq. (JSOj) as 



./() 



e *(-) s -i(l_ s )«-i dS; (35) 



where vE'(s) = 20 (Ins — s). The function \l/(s) has its 
maximum exactly at s = 1, the upper integration limit. 
The largest contribution to the integral comes from the 
small region 0(1/ \/0) near s = 1. Therefore, it suffices 
to expand ^(s) up to the second order in (s— l) 2 , replace 
the factor s _1 by 1 and extend the lower integration limit 
to — oo. The result is 



-2Q 



(1 



s) n - x ds 



e- 2ll T(n/2) 
2 0"/2 



(36) 
(37) 

El. 



Therefore, for n = 0(1), we obtain 

n nl 

which, for n ^> 1, coincides with that given by Eq. 

Figure [2] compares our analytical result (|31[) with (i) 
a numerical solution of the (truncated) master equation 
(P) with (d/dt)P n (t) replaced by zeros and Po = 0, (ii) 
the prediction from the Fokker-Planck equation for this 
problem [Eq. (|B2[) of Appendix B] , and (iii) the gaussian 
distribution ()34|) for [iE\t < 1, In the central part all 



the distributions coincide. The Fokker-Planck approxi- 
mation strongly underpopulates the low-n tail, and over- 
populates the high-n tail. On the contrary, the gaussian 
approximation strongly overpopulates the low-n tail and 
underpopulates the high-n tail. Our analytical solution 
(I3ip is essentially indistinguishable from the numerical 
result, even at small values of n. Actually, it is in good 
agreement with the numerics already for — 0(1), and 
the agreement improves further as O increases. 




FIG. 2: (Color online) The natural logarithms of the ana- 
lytical result (|31[1 for the quasi-stationary distribution (the 
dots), of the distribution obtained by a numerical solution 
of the (truncated) master equation (0 (the solid line), of 
the stationary solution (IB2|) of the Fokker-Planck equation 
(the dashed line), and of the gaussian distribution (|34|l (the 
dashed-dotted line), for Q — 30 and /j,Eit <C 1. 

We also computed the ground-state eigenvalue by solv- 
ing Eq. ([4]) numerically with the boundary conditions 
G(l,t) = 1, d x G(— l,t) = and the initial condition 
G{x,t = 0) = x n °. At times ui > I/O, the numeri- 
cal ground-state eigenvalue Ei um can be found from the 
following expression: 



J3J" 



= -— ln[l 



G num (0,t)} , 



(38) 



where G num (x,t) is the numerical solution for G(x,t), 
and the result in Eq. (|38p should be independent of time. 
A typical example is shown in Fig. d and a good agree- 
ment with the theoretical prediction (|2"Tj) is observed. 



V. STATISTICS OF THE EXTINCTION TIMES 

The quantity Po(t), given by Eq. (|27|) . is the proba- 
bility of extinction at time t. The extinction probability 
density is p(t) = dPo(t)/dt. Using Eq. (|27[) . we obtain 
the exponential distribution of the extinction times: 



p(t) ~ fxEte'^ 11 at At > 1 
The average time to extinction is, therefore, 



tp(t)dt ~ (fiEx)- 1 = 



/iO 3 / 2 



20(l-ln2) 



(39) 



(40) 
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100 



FIG. 3: Shown is the ratio of the numerical ground-state 
eigenvalue E" um from Eq. (|38|) and the approximate ana- 
lytical value of Ei from Eq. (|2T}, for Q, = 20 and no = 100. 
The deviation from 1 is about 5.6 percent, that is within error 
0(1/0). 
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FIG. 4: (Color online) Shown are the extinction probability 
Po(t) [Eq. (|27[) ] (the dashed line) and the numerical solution 
of Eq. (QJ (the solid line) at x = 0, for f2 = 20 and n = 100. 



This is in full a gree ment with the result of Turner and 
Malek-Mansour and in disagreement with the pre- 
diction from the Fokker-Planck approximation, given by 
Eq. (|B5|) . and with the prediction from the gaussian ap- 
proximation, given by Eq. (|B7|) . see Appendix B. 

Figure 2] compares the analytical result (f2"T)l for Po (t) 
with the extinction probability P$ um {t) = G num (0,t) 
found by solving Eq. Q numerically as described at the 
end of the previous section, and a very good agreement 
is observed. 



VI. FINAL COMMENTS 

We calculated, at intermediate and long times, the 
complete probability distribution, including the quasi- 
stationary distribution of the long-lived metastable state, 
and the extinction time statistics in a (non-single-step) 



branching-annihilation reaction. To this end we em- 
ployed the spectral method, recently developed by the 
authors [28l |29||. We also used this example to illus- 
trate how the "lacking" boundary condition of the spec- 
tral method emerges in the theory. 

The spectral method reduces the problem of finding 
the statistics to that of finding the ground-state eigen- 
value and eigenfunction of a linear differential operator 
emerging from the generating function formalism. The 
quasi-stationary distribution that we have calculated an- 
alytically is in excellent agreement with numerics. The 
two widely used "rival" approximations: the Fokker- 
Planck approximation and its reduced version, the gaus- 
sian approximation, perform well only in the peak region 
of the quasi-stationary distribution. They both fail in the 
tails of the distribution and, as a result, cause exponen- 
tially large errors in the estimates of the mean extinction 
time. 

It is worth reiterating that, for single-step birth-death 
systems, the quasi-stationary distribution can be found 
directly from a recursion equation for P n , obtained by 
putting Po — 0, assuming a zero flux into the empty state, 
and replacing (d/dt)P n (t) by zeros in Eq. JT]), see e.g. 
[i~2| . For multi-step systems such recursion equations are 
not generally soluble analytically. 

In conclusion, the spectral method is a powerful tool 
for calculating the quasi-stationary distributions and ex- 
tinction time statistics of a host of multi-step birth-death 
processes possessing a metastable state and an absorbing 
state. 
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Appendix A 

Here we will derive the result given by Eq. (f2"0")) . First, 
we calculate the inner integral 

f°° „ / P \ 2U 

h = J ^(l + r) 2Q e-™ r dr= (^J r(2Q) . (Al) 

As fi ^> 1, we can use the large- argument asymptotics of 
the gamma-function and obtain 



(A2) 



Second, at 3> 1, the integral 



x e 2fls 



o (1 + *) 
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ds 



f e 2nr( - s Us (A3) 
Jo 
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where T(s) = s — ln(l + s), receives the largest contribu- 
tion in the vicinity of s = x (remember that 1 — x <C 1). 
Then, expanding T(s) around s = x, we obtain 

T(s) =x-ln(l + x) + — (s-x) + ... . (A4) 
1 + x 

Extending the lower integration limit to — oo and eval- 
uating the remaining elementary integral we obtain, in 
the leading order, 



h 



i7(i + x) 2n 



(A5) 



Appendix B 

What are the predictions from the Fokker-Planck (FP) 
approximation for the quasi-stationary distribution and 
the mean extinction time of the branching-annihilation 
problem? The FP description introduces an (in general, 
uncontrolled) approximation into the exact master equa- 
tion |T]) by assuming n 3> 1 and treating the discrete 
variable n as a continuum variable. The FP equation 
can be obtained from Eq. ((TJ) by a Kramers-Moyal "sys- 
tem size expansion" P, 0| (in our case, expansion in the 
small parameter 17 _1 <C 1). Using this prescription, we 
obtain after some algebra 

+ ~[2n(2n + n)P(n,t)]} . (Bl) 



receives its main contribution from the vicinity of n = 0. 
Therefore, one can use the saddle point method for the 
inner integral and a Taylor expansion for the outer one. 
The result is 



1 



3 u!7 3 / 2 



■ ln3— l) 



(B5) 



Comparing it with Eq. (|40|) . one can see that the FP ap- 
proximation gives a poor estimate of the mean extinction 
time, as it introduces an exponentially large error. 

The central (gaussian) part of the quasi-stationary dis- 
tribution, |n — f2| <C Q (|B3I) . can be correctly obtained by 
keeping only leading order terms, in the small parameter 
|n — S7|/f2 <C 1, in the FP-equation: 



dP(n,t) u f d 
M ~ 2 \~~dn 

+ i^[6^P(»,«)] 



[2Sl(Sl-n)P(n,t)] 



(B6) 



Indeed, the zero-flux steady state solution of this equa- 
tion yields the gaussian distribution (|B3[) . 

Finally, what would be the prediction for the mean 
extinction time from the reduced FP description, that is 
the one in terms of Eq. (|B6|I ? Here one would obtain 



e3« s dii 



e 3 3f2 

3ui7 2 



-dk 



fin 3 / 



37T n 

e 3 , 



which again gives an exponentially large error as com- 
pared with the accurate result (T4H|) . 



The quasi-stationary distribution of the metastable state 
corresponds to the (zero-flux) steady state solution of the 
FP equation. In the leading order in 1/17 we obtain 



P at (n) ~ (3ttO) 



-1/2 & n-r 



Hi In 



3£2 



(B2) 



where only the central (gaussian) part of the distribution 
contributes to the normalization. In fact, the distribu- 
tion (IB2|1 is accurate only in the peak region, \n— fi| -C 17 
(see Section IV), where it reduces to a gaussian distribu- 
tion with mean 17 and variance 317/2: 



gauss 



(n) = (3tt17) 



-1/2 e -i 



"SIT 



(B3) 



Were the FP equation (|B 1 1) valid for all n, one could use 
it to find the mean time to extinction Tpp (conditional 
on non-extinction prior to reaching the quasi-stationary 
state) by a standard calculation, see e.g. Ref. fl. This 
calculation would yield 



tfp 



2n 



-3H/2 



,-fc 



1 + f) 



^fc(2fc + 0) 



-/." / .. " — — dk. 

(B4) 

As 17 3> 1, the inner integral receives its main contri- 
bution from the vicinity of k = 17. The outer integral 



Appendix C 

Here we calculate the n-th derivative of tpb(x), given 
by Eq. (I16p . at x — 0. The first derivative is 



V ' b (x) = -—1 — 1 



2Q 



■2Q.1 



17 V 217/ {l + x) 2n 
x {r[217 + 1] — T[217 + 1, 217(1 + x)]} . (CI) 

Let us introduce two auxiliary functions: 

{P[217 + 1] - T[2S7 + 1, 217(1 + a;)]} 



h(x) 



(1 + x) 



21 > 



2Q.X 



(C2) 



Using Eqs. (|C1[) and (|C2[) . we can write the n-th deriva- 
tive of tpb(x) [that is, the (n — l)-th derivative of tp' b (x)] 
at x = as 



dx r ' 



_Ei (_^_V il \ - 



x=0 



— (n - 1)! 



17 V217/ 



fc=0 



k\(n - k- 1)! 



x /W(x) ^"-^^(x) , (C3) 

a=0 x=0 
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where f^ k > (x) is the fc-th derivative of f(x), and the same 
notation is used for h(x). 

After some algebra, we find that the fc-th derivative 
(fc > 1) of f(x) at x = is [Mil] 



d k f{x) 



dx k 



-l) fc 2ft [T(2n + fc) - r(2ft + fc, 2ft)] 



(C4) 



(C5) 



Using Eqs. (28]), fC3]) . (|C4|) and (C5]), we obtain for 
n > 1: 



Now, the fc-th derivative of h(x) at a; = is 
d k h(x) 



dx k 



P n (t) = e-^* 



2Ei /_e_\ 2 ° 
~ V2ft/ ^ 



(-l) fe (2ft) r 



-k-l 



k=0 



fcl(n-fc-l)! 



[r(2ft + fc) -r(2ft + fc,2ft)] 



(C6) 



Actually, for n = 1 one has 



P 1 (t) ~ (TT/ft) 1 / 2 



1+0 



(ft- 1 ^) 



and the sub-leading term 0(ft -1 / 2 ) has been neglected 
in Eq. (|C6p . Finally, using Eq. (fT7|) and changing the 
order of summation in Eq. (|C6|) . we obtain the following 
result for n > 1: 



, , 2E 1 (2ft)"- 1 e 2 °r(2ft) , N „ p , 

n 1 (Zil + n) 

(C7) 

where iFi^ajb^x) is the Kummer confluent hypergeomet- 
ric function [35|. To avoid excess of accuracy, we need 
to work with the large-ft asymptotics of this result, see 
Section IV. 
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